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Abstract 

For Ising spin models which bear the spin-ice 
type macroscopic (quasi-)degeneracy, conven- 
tional classical Monte Carlo (MC) simulation 
using single spin flips suffers from dynamical 
freezing at low temperatures (T). A similar 
difficulty is seen also in a family of Heisen- 
berg spin models with easy-axis anisotropy or 
biquadratic interactions. In the Ising case, 
the difficulty is avoided by introducing a non- 
local update based on the loop algorithm. We 
present an extension of the loop algorithm to 
the Heisenberg case. As an example of its 
application, we review our recent study on 
spin-glass (SG) transition in a bond-disordered 
Heisenberg antiferromagnet coupled to local 
lattice distortions. 

1 Introduction 

Recently, increasing attention has been de- 
voted to low-temperature behavior of geomet- 
rically frustrated magnets [1]. Spin glass, in 
which spins are frozen randomly, is one of low- 
T phases widely observed in geometrically frus- 
trated materials. However, it is unclear so 
far how the nature of SG is different from the 
canonical one driven solely by randomness. 

An antiferromagnet on a pyrochlore lattice 
(Fig. d]) is a typical example of geometri- 
cally frustrated spin systems. Recently, sev- 
eral puzzling SG behaviors have been pointed 
out experimentally in pyrochlore-based mag- 
nets. One of the surprising aspects is that. 



in these SG materials, the SG transition tem- 
perature Tf appears to be almost indepen- 
dent of the strength of disorder A; e.g., 
for (La^Yi_^)2Mo207, Tf - 22K stays al- 
most constant for x < 0.5 [2J. Similar 
plateau behavior of Tf is also observed in 
(Zni_a;Cda;)Cd204 [3J. Another distinctive as- 
pect is that, for these SG materials, Tf is much 
higher than a theoretically expected value for 
a moderate strength of disorder A [U E]. 
These behaviors suggest that some important 
factor is missing in the previous SG theories: 
A candidate is the magnetoelastic coupling. 
For example, various microscopic probes have 
pointed out importance of local lattice distor- 
tions in Y2M02O7 although this material ex- 
hibits no uniform lattice distortion. They are 
crucial also in (Zni_a;Cda;)Cr204 because it ex- 
hibits the spin-lattice ordering at x = 0. 

Motivated by these puzzles, we recently in- 
vestigated effects of the magnetoelastic cou- 
pling on the spin-glass transition by consid- 
ering the following classical Heisenberg spin 
model [7J: 



rt — y ^ Jij^i ' ^j ~ bij ybi 



Sjf 



(1) 



where Si denotes a Heisenberg spin at site i 
and the sum runs over nearest-neighbor bonds. 
The biquadratic interaction bij (= 6J^j, b > 0) 
is induced by the spin-lattice coupling to local 
lattice distortions [8J. Note that such 'ferro'- 
type biquadratic interaction favors collinear 
spin configurations. We introduce static bond 



disorder as a uniformly-distributed random- 
ness as Jij G [J- A, J+A] with < A < J. At 
A = 0, this model exhibits a nematic transi- 
tion at Tc ^ 6, below which spins select a com- 
mon axis without selecting their directions on 
it; the system remains magnetically disordered 
down to zero T. The ground-state degener- 
acy is equivalent to that of a nearest-neighbor 
Ising antiferromagnet on a pyrochlore lattice 
[see Fig. [U^a)]. The ground-state degenerate 
manifold of this Ising model is identified by 
a set of local constraints enforcing two spins 
pointing up and two spins pointing down in 
every tetrahedron [9j. This is called the ice 
rule because of an analogy to the constraint on 
positions of protons in hexagonal ice [TUl ITT] . 
Similar situation is seen also in the so-called 
spin-ice model in which Ising spins along the 
local (111) axes interact with each other ferro- 
magnetically [121 [13] [see Fig. [Tj^b)]. In these 
systems, the degenerate configurations are sep- 
arated by energy barriers of the order of the 
dominant interaction scale J, and the standard 
single-spin-flip algorithm does not work at low 
T ^ J in classical MC simulation. The dif- 
ficulty remains even when the Ising discrete- 
ness is relaxed and spins can fluctuate as in 
the present case of the biquadratic interaction, 
as long as the ground-state manifold retains a 
multivalley structure. Indeed, the single-spin- 
flip algorithm does not work at low T ^ 6 for 
the model ([T]) even at A = 0. A similar situa- 
tion is seen also in classical Heisenberg models 
with single-ion easy-axis anisotropy. 

In the Ising case, the difficulty is avoided by 
introducing a global flip called the loop flip, in 
which one reverses all Ising spins on a specific 
closed loop passing through tetrahedra [14j; 
the loop is chosen so that the spins are up 
and down (or inward and outward) alterna- 
tively along the loop as illustrated in Fig. [H 
The loop flip connect different ice-rule states 
bypassing the energy barriers. In the Heisen- 
berg case, however, it is nontrivial how to de- 
fine the loop with alternating spins. Moreover, 
the loop flip procedure is not unique because 
of the continuous degrees of freedom. These 
argued us to extend the loop algorithm to the 
Heisenberg case to investigate low-T properties 
of the model ([T]). 




Figure 1: A 16-site cubic unit cell of the 
pyrochlore lattice is shown with Ising spins 
along a global axis [(a)] and local (111) axes 
[(b)]. Spins are denoted by arrows, (a) White 
circles represent spins pointing upward, while 
black circles the opposite. The spin configura- 
tion is an example of the spin-ice type states. 
The hexagon with a bold dashed line denotes 
an example of loops with alternating black 
and white sites, (b) The ice-rule configura- 
tion equivalent to (a) is shown with Ising spins 
along the local (111) axes. A and B represent 
two different types of tetrahedra. Black cir- 
cles represent spins pointing inward in terms 
of type-A tetrahedra, while white circles the 
opposite. 

This report is organized as follows. In § O 
we briefly review the loop algorithm for Ising 
models [H]. Section [3] is devoted to the exten- 
sion of the loop algorithm to Heisenberg spin 
systems and its benchmarks [151 [IS]- In § [H 
we report results of our numerical study on 
the model ([T]) using the extended algorithm. 

2 Loop algorithm for Ising 
spin systems 

Before considering an extension of the loop al- 
gorithm to Heisenberg spin systems, here we 
briefly review the loop algorithm for Ising spin 
models [1]. To generalize the following dis- 
cussion, we assign black and white to two de- 
grees of freedom of Ising spins in an appropri- 
ate manner. For example, black and white sim- 
ply correspond to up and down spins, respec- 
tively, for the antiferromagnetic Ising model 
[Fig. [IJa)]. For the spin ice model [Fig. Hb)], 
black and white represent inward and outward 
spins in terms of type-A tetrahedra, respec- 
tively. Then the loop flip consists of two steps; 



first, we identify a closed loop which consists of 
alternating alignment of black and white sites, 
and then we try to flip all Ising spins on the 
loop. 

Such a closed loop can be constructed by us- 
ing the short loop algorithm |14| . In the short 
loop algorithm, one traces a path through 
alternating black and white sites in ice-rule 
tetrahedra. A loop is formed when the path 
encounters any tetrahedron already included 
in the path as illustrated in Fig. O At finite 
T, thermal fluctuations induce "defect tetra- 
hedra" in which the ice-rule condition is vio- 
lated. To maintain detailed balance, the path 
must be traced so that it does not involve de- 
fect tetrahedra on it. 

After the construction of a closed loop, all 
colors on the loop are reversed simultaneously 
by flipping the spins. When all the ice-rule 
states are energetically degenerate, the flip is 
always accepted (rejection free) in the MC 
samphng because the flip does not the total 
energy. When there are residual interactions 
which lift the degeneracy, the loop flip is ac- 
cepted according to the Metropolis criterion. 

At finite T, the loop flip update does not sat- 
isfy ergodicity because it changes neither the 
spin configurations in defect tetrahedra nor the 
number of defect tetrahedra. It is, therefore, 
necessary to use the loop flip together with an- 
other update such as the standard single-spin 
flip for retaining the ergodicity. 



Defect tetrahedron 




Figure 2: Schematic picture for a loop con- 
struction by tracing a path through ice-rule 
tetrahedra. The path is made of alternating 
black and white sites. For simplicity, the fig- 
ure shows a (111) kagome layer with connected 
tetrahedra. The path is denoted by a dashed 
line, and its left part represents an example of 
a closed loop. 



3 Extended Loop algorithm 

Now we extend the loop algorithm to Heisen- 
berg spin systems with the spin-ice type de- 
generacy: Heisenberg spin systems with (1) 
single-ion anisotropy [E] and (2) biquadratic 
interactions [IH]. Section ED is devoted to an 
overview of the extended loop algorithm. In 
§ 13.21 and § 13.31 we review the detailed proce- 
dure of the algorithm for the cases (1) and (2), 
respectively. Benchmark results are also given. 

3.1 Overview of the extended algo- 
rithm 

The extended loop algorithm consists of the 
following three steps: 

1. We first project the Heisenberg spin Si 
onto an appropriate projection axis di to 
assign black and white colors at every site 
I. 

2. Then we construct a loop consisting of al- 
ternating black and white sites. 

3. All colors on the constructed loop are re- 
versed simultaneously. 

In the case of the single-ion anisotropy, the 
projection axis di is simply given by the easy 
axis at site i. For example, in the case of an- 
tiferromagnets with easy-axis anisotropy along 
the z axis, which are extensions of the Ising an- 
tiferromagnet [Fig. [T]^a)], we set di — (0,0,1) 
for all the sites. While, in the case of fer- 
romagnets with the local (111) anisotropy, 
which are natural extensions of the spin ice 
model [Fig. U^b)], we set di to the direction 
connecting the centers of neighboring tetrahe- 
dra from type B to A. In the case of the bi- 
quadratic interactions, however, the systems 
retain 0(3) spin rotational symmetry and have 
no explicit anisotropy axis to project the spins 
on. Therefore, it is necessary to deduce the 
common axis selected by spins for each MC 
sample. In § 13. 3[ we introduce a simple way to 
determine the projection axis. 

Once the projection axis di is defined at ev- 
ery site, we assign black and white colors to 
sites at which S^ • 5^ > and Si - di < 0, re- 
spectively. Based on this definition, we can 



construct a closed loop with alternating black 
and white sites by following the short loop al- 
gorithm similarly to the Ising case. 

In the step 3, as mentioned above, the way 
to reverse black and white is not unique be- 
cause of the continuous degrees of freedom. 
Three different ways are illustrated in Fig. [3l 
(1) flip xyz^ (2) flip parallel and (3) rotate. In 
flip xyz, all three Cartesian components of Si 
are reversed as S^ ^ — S^, while in flip parallel^ 
only components parallel to the easy axes, S'^y, 
are reversed as Si ^ Si — 2(5^ • 5^)5^. In ro- 
tate, which is applicable to systems with global 
anisotropy axes, one translates every spin to 
the neighboring site on the loop simultaneously 
in the same direction. 

For models which retain the ground-state 
spin-ice type degeneracy, one might expect 
that these updates become equivalent and al- 
ways accepted at low T; this is naively ex- 
pected since thermal fluctuations vanish and 
all the ice-rule conflgurations with spins par- 
allel to the easy axes become energetically de- 
generate. However, this is not the case: As 
discussed in the following sections, careful con- 
sideration on the energy change is necessary to 
choose an eflicient method. 



(2) flip paralleU;.f 




(3) rotate 
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Figure 3: Different ways to reverse black and 
white: (1) flip xyz, (2) flip parallel, and (3) 
rotate, di is the projection axis at site i. 
The rotate is applicable to systems with global 
anisotropy axes. See the text for details. 

Now we comment on technical aspects of 
implementation of the algorithm. We imple- 



mented a computation code based on the ex- 
tended loop algorithm, single-spin-flip algo- 
rithm and exchange MC method [I7|, which 
is used in the following MC simulations. As 
illustrated in Fig. [H one MC step consists of 
a sweep of the lattice by sequential single-spin 
flips, followed by the loop update and replica 
exchange between neighboring temperatures. 
The loop flips are repeated until the number 
of tetrahedra visited in the loop construction 
exceeds the number of lattice sites. The loop- 
flip section takes CPU time comparable to the 
single-spin-flip sweep. 






o 



Sequential single-spin flips 



Loop update 

0. Estimation of the common axis 
(only for biquadratic interactions) 

1 . Color assignment (black/white) 

2. Loop construction ■ 

3. Loop flip using flip xyz, flip parallel, 
or rotate 



Replica exchange, sampling of physical 
quantities, etc. 



Figure 4: Flowchart of our computational code 
with the extended loop update, the single-spin 
update, and the replica exchange MC method. 



3.2 Heisenberg spin systems with 
single-ion anisotropy 

In this section, we review the loop algorithm 
extended to Heisenberg models with easy-axis 
anisotropy p^. We start with a simple Hamil- 
tonian: 



n = jJ2Si-s,-D,J2 



Oii 



(2) 



{hj) 



where Si denotes a classical Heisenberg spin at 
site i (we take l^^l = 1) and Di (> 0) is the 
single-ion easy-axis anisotropy. The easy axis 



Si (\ai\ = 1) defines the projection axis in the 
loop update. Although the exchange interac- 
tion is limited to nearest neighbors for simplic- 
ity here, the following algorithm is applicable 
to more general models with farther-neighbor 
or bond-dependent interactions. We consider 
periodic systems of cubic geometry with L^ 
unit cells with totally Ng = 16L^ spins. We 
take the energy unit as | J| = 1. 

For models which retain the spin-ice type 
degeneracy in the ground state, spins fluctu- 
ate around the easy axes by angles of 0{-\/T) 
at low T <^ Di. Under the influence of the 
thermal fluctuations, the energy change in the 
loop flip is estimated as followq^: 



AEocT 
AEocT^ 



{flip xyz), 
{flip parallel) . 



Because the acceptance rate of the Metropolis 
algorithm is given by min{l,exp(— A£'/T)}, 
thermal fluctuations are irrelevant 
for flip parallel in the sense that 
limT^oexp(— A£^/T) = 1. On the con- 
trary, thermal fluctuations are relevant for 
flip xyz since lim^^^o exp(— AE'/T) < 1. 
This consideration indicates that flip parallel 
becomes refection free in the limit of T ^ 0, 
while flip xyz not even for models which retain 
the spin-ice type degeneracy. 

The efficiency of the loop flip is demon- 
strated in Fig. [5]^a) for the model ([2j) with the 
antiferromagnetic exchange interaction J = 1 
and Si = (0, 0, 1). This model retains the spin- 
ice type degeneracy in the ground state. At 
low T <^ I J| = 1, spin configurations are grad- 
ually enforced to satisfy the ice rule, and the 
acceptance rate of the single-spin flip, Pgingie^ is 
suppressed below T ^ \J\ and vanishes in the 
low-T limit. On the contrary, the acceptance 
rate of loop flips increases at low T. As shown 
in Fig. [5]^a), the probability that a closed loop 
is successfully formed, Pioop? steeply increases 
below T r^ I J|, indicating that almost all tetra- 
hedra start to follow the ice rule below this 
temperature. At the same time, the accep- 
tance rate of flips of a formed loop gradually 
increases at low T < \J\ and remains finite; 



here, Pxyz and Pparaiiei are the rate for flip xyz 
and flip parallel^ respectively. The total ac- 
ceptance rate of the loop flip is given by the 
product as Pioop x P^yz or Pioop x Pparaiiei, and 
it sharply increases at T < |J|, compensating 
the decrease of Pgingie- 

As clearly indicated in Fig. [5]^a), the accep- 
tance rate of flip parallel is always larger than 
that oi flip xyz, being consistent with the above 
argument. In particular, Pparaiiei approaches 1 
(rejection free) as T ^ 0, whereas Pxyz goes to 
a smaller value ^0.5. The reduction of Pxyz 
becomes larger for smaller anisotropy Di. This 
is demonstrated at T = 0.1 in Fig. [5]^b); Pxyz 
decreases almost exponentially with 1/Di. On 
the other hand, Pparaiiei is almost independent 
of Di and remains rejection free at T ^ in 
the wide range of Di. 

3.3 Heisenberg spin systems with bi- 
quadratic interaction 

In this section, we review the loop algorithm 
extended to classical antiferromagnetic Heisen- 
berg models with biquadratic interactions. We 
start with a Hamiltonian of a simple form: 



n 






Si • Sn] -b{ Si- S. 



-. ^ \ 2 



(3) 



^For more detailed discussion, refer to ref. [15] and 
a short note in reference 22 of Ref. [161 • 



where 6 (> 0) is the biquadratic interaction. 
The model ([T]) reduces to this model when 
A = 0. We note that such 'ferro'-type bi- 
quadratic interaction originates in quantum 
and thermal fluctuations as well as the spin- 
lattice coupling. We consider the antiferro- 
magnetic exchange interaction J > 0, and take 
the energy unit as J = 1. The sum runs 
over nearest-neighbor bonds. The following al- 
gorithm is applicable to more general models 
with fart her- neighbor or bond-dependent in- 
teractions such as the model ^ with A > 0. 

As mentioned above, it is necessary to de- 
duce the common axis Q selected by b for each 
MC sample. Here, we explain a simple way to 
determine the projection axis. We first pick up 
a set of A^T tetrahedra {Tm} {m — 1,- • • , A^t) 
randomly from the whole system. Starting 
from an initial guess Sq [we take Sq = (0, 0, 1)], 
the normalized projection axis a is obtained it- 



(a) 1.0 




Figure 5: (a) Temperature dependences of the 
acceptance rates of the single-spin flip (Psingie)^ 
the probabihty of formation of closed loops 
(Pioop), the acceptance rates of flip of a formed 
loop by flip xyz {Pxyz) and by flip parallel 
(Pparaiiei)- The data are calculated for the 
model ([2]) at i^i = 5.0 with the antiferromag- 
netic exchange interaction J — 1 and di — 
(0,0, 1). The data for L = 2 and L = 4 are 
denoted by crosses and fllled squares, respec- 
tively, (b) Di dependence of Pparallel, Pxyz, 

Pioop, and Psingie at T = 0.1 for L = 2. 



eratively by 



«n+l O^ 



E 



s\gn{Si • dn)Si. 



Here the sum is taken over all spins belong- 
ing to the selected tetrahedra {7^}, and n {— 
0, 1,--- ,nmax — 1) is the index of the itera- 
tion. For larger N^ and rimax, the resultant 
a — Srimax gives a better approximation of Q. 
In practice, we take A^t = 16 and rimax = 6 
for the system sizes L > 3 in the following MC 
simulations. It should be noted that, to en- 
sure the detailed balance, loops must be con- 
structed avoiding the tetrahedra included in 
{Tm} as well as defect tetrahedra in which the 



ice rule is violated: Otherwise, the loop flip 
becomes irreversible because the flip changes 
a. Because loop flips do not change 5, a is 
determined once at the beginning of each MC 
step as shown in Fig. [H The computational 
cost for estimating a is negligible in practical 
calculations. 

Now we show benchmark results in Fig [6l 
The number of spins in the system Nq is given 
by 16L^. For 6 > 0, the model ([3D exhibits 
a nematic transition at Tc ^ 6, below which 
spins select a common axis. At low T com- 
pared to h and J, spin conflgurations are en- 
forced to satisfy the 'two-up two-down' ice rule, 
and the acceptance rate of the single-spin flip, 
^single 5 is Suppressed. This is demonstrated in 
Fig. EKa) with h = 0.2. While, the probability 
that a closed loop is successfully formed, Pioop? 
steeply increases below Tc ^ 6, indicating that 
almost all tetrahedra start to follow the ice rule 
below Tc. The acceptance rate of flips of a 
formed loop also increases below T^ and re- 
mains flnite as T ^ 0; here, Pxyz, ^parallel? and 
^rotate are the rates for flip xyz^ flip parallel, 
and rotate, respectively. The acceptance rate 
of the loop flip sharply increases at T < Tc, 
compensating the decrease of Psingie- 

In contrast to the case of the single-ion 
anisotropy, the most eflicient loop flip depends 
on the value of b as demonstrated in Fig. [6](b): 
Pparallel bccomes most eflicient as 6 ^ 0, while, 
in the opposite limit (i.e., b -^ +oc), flip xyz 
becomes rejection free but the other two not. 
In the intermediate regime, i.e, 0.1 < 6 < 0.5, 
rotate is most eflicient. The difference of the 
efficiency of the loop ffips are understood by 
the following consideration. Considering a 
given state at a finite T well below Tc, its en- 
ergy measured from the ground-state energy 
is given by £^ = £^j + £^5 = 0(T), where Ej 
and El, are the energies corresponding to the 
first and second terms in eq. (Ej), respectively. 
The three loop fiips change the two contribu- 
tions in different ways. The flip xyz conserves 
£■5, while the other two not. This is why the 
flip xyz becomes most efficient at 6 ^ oc where 
E '2:^ El) and \Ei)\^ |^j|- I^i the opposite limit, 
the contribution of Ej becomes dominating in 
the energy change because \Ej\ ^ \Ei)\- The 
flip xyz, flip parallel and rotate change Ej by 



0(r), 0(r2) and 0{T), respectively. This is 
why flip parallel is most efficient at smah b. 
The complete argument is given in ref. [16j. 
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Figure 6: (a) Temperature dependences of the 
acceptance rates. For the definition of the ac- 
ceptance rates, see the text. The data are 
calculated at 6 = 0.2 and L = 8. The verti- 
cal broken line denotes the nematic transition 
temperature Tc. (b) b dependence of the ac- 
ceptance rates at T = 0.02. The most efficient 
method depends on the value of b. 



4 Application 

Now, we review results of our MC study on 
the model ^ utilizing the extended loop al- 
gorithm and discuss the obtained A-T phase 
diagram. When 6 = 0, the model reduces to 
the previously studied one [H El El [IE] • In this 
case, by turning on A, the SG transition ap- 
pears roughly linearly as Tf 2^ O.IA [HEIEI. 
As already mentioned, for a finite spin-lattice 
coupling 6 > 0, the present model exhibits a 
nematic transition at Tc ^ 6 for A = 0. Our 
interest is how the SG transition appears in 
the presence of b by turning on A. 

In the following MC simulations, we take 
6 = 0.2 and employ rotate^ which is the most 
efficient flip at 6 = 0.2. We also adopt the ex- 
change MC method ^^ and the overrelaxation 



update [19J to further accelerate MC dynam- 
ics. To identify the SG transition, we calculate 
the SG susceptibility XSG = ^s^eA' where the 
Edwards- Anderson order parameter g^^ [20] 
is defined as the overlap of two independent 
replicas with the same interaction set {Jij}. 

The computational code is parallelized by 
using the message passing interface (MPI) li- 
brary. The MPI processes are divided into in- 
dependent groups which handle different {Jij}. 
In the replica exchange, spin configurations are 
swapped within a process or across neighbor- 
ing processes. In the following simulation, we 
take four independent replica for each interac- 
tion set. For L = 5, a MC run of 1.2 x 10^ steps 
with one interaction set and 64 temperatures 
takes 10 fours using 64 CPU cores in System B 
of the Supercomputer Center (ISSP, Univ. of 
Tokyo). 
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Figure 7: (a) Comparison of thermalization 
processes of g^^ with/without the loop update 
{rotate), (b) Calculated A-T phase diagram of 
the model (|I]). 

In Fig.[71^a), we compare thermalization pro- 
cesses of q^js^ with/without the loop update 
at A = 0.1 and T = 0.08 for the system 
size L — 2 with 128 spins. For the both 
cases, we took 16 temperatures in the range 



of 0.08 < T < 0.2 for parallel tempering. 
The MC dynamics without the loop update 
exhibits a severe freezing, and does not reach 
thermal equilibrium after 2 x 10^ MC steps. 
However, the thermalization process is greatly 
accelerated by the loop update and g^^ con- 
verges to the thermal-equilibrium value after 
2.5 X 10^ MC steps. These clearly show the 
advantage of the loop algorithm in investigat- 
ing the low-T properties of the present model. 
The advantage becomes more pronounced for 
larger L, because the number of the ice-rule 
states grows as L increases. 

Figure [71(b) presents the phase diagram ob- 
tained by MC simulation with the loop up- 
date for L < 5. By introducing the disorder 
A, the SG transition appears at a finite Tf. 
In the weakly-disordered region, i.e., A < 6, 
Tf is roughly proportional to A as Tf 2:^ A. 
A remarkable point is that Tf is largely en- 
hanced by 6 [1[ El El : The enhancement factor 
is, e.g., about 5-10. At A 2^ 6, Tf appears 
to merge into T^ with showing multicritical 
behavior. For larger A, Tf (= Tc) becomes 
nearly independent of A, being in sharp con- 
trast to the previously-reported SG behavior, 
i.e., Tf oc A m El |6]. These peculiar SG be- 
haviors are ascribed to reduced thermal fluc- 
tuations in the semi-discrete degenerate mani- 
fold emergent below T^. The plateau behav- 
ior of Tf at a largely enhanced value gives 
an explanation for the puzzling behaviors in 
the pyrochlore-based antiferromagnets such as 
(La^Yi_^)Mo207 and (Zni_^Cd^)Cd204. 



5 Summary 

We have reviewed an extension of the loop al- 
gorithm to Heisenberg spin systems. In § El 
we have explained detailed procedure of the ex- 
tended loop algorithm and demonstrated its ef- 
ficiency for Heisenberg spin models with single- 
ion anisotropy and those with biquadratic in- 
teractions. Finally, we have reported results 
of our recent numerical study on spin-glass 
transition in a bond-disordered antiferromag- 
net coupled with local lattice distortions. 

This report is based on the recent works 
done with Y. Tomita and Y. Motome. 
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